function FVSWE2D 
% File name: FVSWE2D.m 
% This m-file is associated with the FREE textbook:  
%             'Introductory Finite Volume Methods for PDEs'  
% which may be downloaded from bookboon.com. 
% 
% Notes:  
%  
% 1. This code uses the 2D Lax-Friedrichs solver and is an extension 
% of previous 1D codes by Clive Mingham. 
% 
% 2. Results are compared to the exact solution of the 1D dam break 
% (Riemann) problem. 
% 
% Description:  Solves the system of 2D PDEs, dU/dt + div(H) = 0, 
% using a finite volume scheme on a 2D Cartesian mesh.   
% 
% This equation corresponds to the finite volume form: 
% 
% doubleintegral(dU/dt dA) + lineintegral(H.n ds) = 0. 
% 
% In this case the equations are the 2D shallow water equations (SWE) so:  
% Water depth: h=h(t,x,y) [m]. 
% Water velocity: v = vx i + vy j where vx=vx(t,x,y), vy=vy(t,x,y) [m/s]. 
% Acceleration due to gravity: g=9.81 [m/s^2]. 
% Geopotential: phi = gh.  
% Matrix of conserved variables: U = [  phi 
%                                     phi vx 
%                                     phi vy ]. 
% 
% Flux density: H =[    phi*vx       i  + [     phi*vy         j 
%                   phi*vx^2+phi^2/2           phi*vx*vy 
%                      phi*vx*vy   ]        phi*vy^2+phi^2/2 ] 
% 
% subfunctions: freadmesh, fcellareas, fsidevectors, finitialU,  
%               finterfaceflux, fcalcdt, fsolver,fplotresults, fgetHfromU, 
%               fcelldistances, fdisplay, fdrawmesh. 
% 
% Verification:  Results look similar to the exact solution.  The code  
%                really ought to be verified by a pen and paper calc. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% 
clc; clf; clear all; % clear everything 
g=9.81; % acceleration due to gravity [m/s^2] 
runtime=1.2;  % runtime [s] 
t=0;          % current time [s] 
timelevel=0;  % current time level 
[x,y,xcen,ycen]=freadmesh; % gets coordinates of the computational mesh 
                           % and cell centres.  Each array contains left 
                           % and right ghost cells so that computational  
                           % indices start at i = 2 and j = 2 and end at 
                           % i = NI+1 and j= NJ+1 for xcen and ycen. 
                           % 
% Dinterfaces holds distances from cell centres to the 4 interfaces 
% Dcentres holds distances from cell centres 4 neighbourng centres 
[Dinterfaces,Dcentres]=fcelldistances(x,y,xcen,ycen); 
%  fdrawmesh(x,y) % for checking small mesh visually 
A=fcellareas(x,y);   % compute and store areas for each cell.  Includes 
                     % ghost cells around the computational domain. 
%area=fdisplay(A) 
S=fsidevectors(x,y); % compute and store cell side vectors in an 
                     % (NI+2)x(NI+2)x4x2 array which includes ghost cells  
                     %  around the computational domain. 
  
% 
[U,hL,hR,vxL,vxR]=finitialU(xcen,ycen); % Put initial cell centre values of  
                        % U in a 3x(NI+2)x(NJ+2) array with ghost cells and 
                        % also stores left and right states for dam break. 
% 
Unext=zeros(size(U)); % Creates an array for the updated U values  
[d,M,N]=size(U);        % d is the number of PDEs (=rows of U) 
d 
NI=M-2 % number of computational cells in the x direction 
NJ=N-2 % number of computational cells in the x direction 
%% 
disp('time marching starts')  
while(t<runtime)  
 timelevel=timelevel+1  
 U=fbcs(U);% Insert boundary conditions into ghost cells at each end 
           % of the domain. U becomes 3x(NI+2)x(NJ+2) so each computational  
           % cell is indexed i=2 to NI+1, j=2 to NJ+2. 
% 
 dt=fcalcdt(A,S,U); % Finds stable time step before each iteration. 
 t=t+dt;  % updates time   
 if (t>runtime) % ensures program stops at exactly t=runtime [s] 
     dt=dt-(t-runtime); 
     t=runtime; 
 end 
 t   % Outputs current time to the Command Window 
 dt 
 Unext=fsolver(U,A,dt,S,Dinterfaces,Dcentres); % implements a time step 
% 
 U=zeros(d,M,N); % create correct sized array of zeros 
 U(:,2:NI+1,2:NJ+1)=Unext(:,2:NI+1,2:NJ+1); % update computational U values 
                                            % before next time iteration 
end  % of while loop 
disp('2D Lax-Friedrichs finite volume scheme ended at time t') 
t 
U=U(:,2:NI+1,2:NJ+1); % remove ghost cells 
xcen=xcen(2:NI+1,2:NJ+1); % remove ghost cells 
ycen=ycen(2:NI+1,2:NJ+1); % remove ghost cells 
%fplotresults(xcen,ycen,U) % plot results 
disp('see plots for comparisons with exact solutions') 
exacth=zeros(NI,1); % exact solution array for water depth 
exactvx=zeros(NI,1); % exact solution array for water speed 
[exacth,exactvx]=fexactSWE1D(runtime,xcen,hL,hR,vxL,vxR); 
subplot(2,1,1),plot(xcen,exacth,xcen,U(1,:,3)/g,'--') 
title('1D dam break: exact (solid) and numerical (dashed) water depths') 
xlabel('x [m]') 
ylabel('depth [m]') 
subplot(2,1,2),plot(xcen,exactvx,xcen,U(2,:,3)./U(1,:,3),'--') 
title('1D dam break: exact (solid) and numerical (dashed) water speeds') 
xlabel('x [m]') 
ylabel('speed [m/s]') 
disp('finished main program, see figures') 
end % FVSWE2D 
%%%%%%%%%%%%%%%%%%%%%%%%%%% subfunctions  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [x,y,xcen,ycen]=freadmesh 
% Verification:  Verified. 
% The computational mesh has NI cells in the x direction and NJ cells  
% in the j direction. 
% The x and y coordinates of the lower left hand corner of cell (i,j) are 
% held in arrays x and y respectively which are both (NI+1)x(NJ+1).  In 
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)), 
% (x(i+1),y(j+1)) and (x(i),y(j+1)).   
% The mesh starts at the origin so x(1,1)=0=y(1,1). 
% Ghost cell vertices are created by reflection of neighbouring 
% cell veritices. 
%  
NI=100;  % number of cells in the i direction 
NJ=5;  % number of cells in the j direction 
x=zeros(NI+1,NJ+1); % allocate correct sized array for cell vertices  
                    % initially without ghost cells. 
y=zeros(NI+1,NJ+1); % allocate correct sized array for cell vertices 
                    % initially without ghost cells. 
% 
dx=1.0; % cell length in x direction 
dy=2.0; % cell length in y direction 
% create uniform mesh 
for i=1:NI+1 
    for j=1:NJ+1 
        x(i,j)=(i-1)*dx; 
        y(i,j)=(j-1)*dy; 
    end % of j loop 
end % of i loop 
% Create ghost cells around whole mesh. 
% First embed x and y in larger arrays 
tempx=zeros(NI+3,NJ+3); 
tempx(2:NI+2,2:NJ+2)=x; 
x=tempx; 
clear tempx 
tempy=zeros(NI+3,NJ+3); 
tempy(2:NI+2,2:NJ+2)=y; 
y=tempy; 
clear tempy 
% Now fill in ghost cell vertices by reflecting cell vertices 
% i.e. if z is the vector from vertex (i,2) to vertex (i,3) then  
% ghost vertex (i,1) has position vector (i,2)-z = 2(i,2)-(i,3) 
for i=2:NI+2 
    x(i,1)=2*x(i,2)-x(i,3); 
    y(i,1)=2*y(i,2)-y(i,3); 
    x(i,NJ+3)=2*x(i,NJ+2)-x(i,NJ+1); 
    y(i,NJ+3)=2*y(i,NJ+2)-y(i,NJ+1); 
end 
% 
for j=2:NJ+2 
    x(1,j)=2*x(2,j)-x(3,j); 
    y(1,j)=2*y(2,j)-y(3,j); 
    x(NI+3,j)=2*x(NI+2,j)-x(NI+1,j); 
    y(NI+3,j)=2*y(NI+2,j)-y(NI+1,j); 
end 
% Fill in the 4 corners by reflecting the diagonals 
% Lower left corner 
v1=[x(2,2),y(2,2)]; % position vector of computational corner 
v2=[x(3,3),y(3,3)]; % position vector of point on interior diagonal 
d=v2-v1; % diagonal vector 
vc=v1-d; % postion vector of reflected ghost corner 
x(1,1)=vc(1); 
y(1,1)=vc(2); 
% Upper right corner 
v1=[x(NI+2,NJ+2),y(NI+2,NJ+2)]; % position vector of computational corner 
v2=[x(NI+1,NJ+1),y(NI+1,NJ+1)]; % position vector of point on interior diagonal 
d=v2-v1; % diagonal vector 
vc=v1-d; % postion vector of reflected ghost corner 
x(NI+3,NJ+3)=vc(1); 
y(NI+3,NJ+3)=vc(2); 
% Lower right corner 
v1=[x(NI+2,2),y(NI+2,2)]; % position vector of computational corner 
v2=[x(NI+1,3),y(NI+1,3)]; % position vector of point on interior diagonal 
d=v2-v1; % diagonal vector 
vc=v1-d; % postion vector of reflected ghost corner 
x(NI+3,1)=vc(1); 
y(NI+3,1)=vc(2); 
% Upper left corner 
v1=[x(2,NJ+2),y(2,NJ+2)]; % position vector of computational corner 
v2=[x(3,NJ+1),y(3,NJ+1)]; % position vector of point on interior diagonal 
d=v2-v1; % diagonal vector 
vc=v1-d; % postion vector of reflected ghost corner 
x(1,NJ+3)=vc(1); 
y(1,NJ+3)=vc(2); 
% 
xcen=zeros(NI+2,NJ+2); % allocate correct sized array for cell centres 
ycen=zeros(NI+2,NJ+2); % allocate correct sized array for cell centres 
% find cell centres by averaging coordinates of surrounding vertices  
% (this works for a general structured mesh and includes ghost cells) 
for i=1:NI+2 
    for j=1:NJ+2 
     xcen(i,j)=(x(i,j)+x(i+1,j)+x(i+1,j+1)+x(i,j+1))/4; 
     ycen(i,j)=(y(i,j)+y(i+1,j)+y(i+1,j+1)+y(i,j+1))/4; 
    end % of j loop 
end % of i loop 
end % freadmesh 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function A=fcellareas(x,y) 
% Verification: Verified. 
% General way to find quadrilateral cell areas. 
% Defines A to include left and right ghost cells. 
 [m,n]=size(x); 
 NI=m-3; % number of computational cells in the i direction 
 NJ=n-3; % number of computational cells in the j direction 
 A=zeros(NI+2,NJ+2);  % allocate correct size array including ghost cells 
                      % around whole computational domain. 
 for i=1:NI+2 
     for j=1:NJ+2 
         % write the diagonal vectors as 3D vectors for the cross product 
         d1=[x(i+1,j+1)-x(i,j),y(i+1,j+1)-y(i,j),0];  
         d2=[x(i,j+1)-x(i+1,j),y(i,j+1)-y(i+1,j),0]; 
         d=cross(d1,d2);    
         A(i,j)=0.5*norm(d); 
     end % of j loop 
 end % of i loop 
end % of fcellareas 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function S=fsidevectors(x,y) 
% Verification: Verified. 
% 
% Finds the 4 side (S) vectors for each cell.   
% Cell sides are numbered 1 to 4 anti-clockwise.  For cell (i,j) side 1  
% is the vector s1 from (x(i,j), y(i,j)) to (x(i+1,j),y(i+1,j)). 
% if s1=<a,b> then the corresponding side vector S1=<b,-a>. 
% The side vector array contains ghost cells. 
% 
% Note that vectors are 1 by 2 arrays so care must be taken because side 
% vectors are stored in higher dimensional arrays. 
% 
[m,n]=size(x); 
NI=m-3; % number of computational cells in the i direction 
NJ=n-3; % number of computational cells in the j direction 
S=zeros(NI+2,NJ+2,4,2); 
for i=2:NI+1 
    for j=2:NJ+1 
        s1=[x(i+1,j)-x(i,j),y(i+1,j)-y(i,j)]; 
        S(i,j,1,1)=s1(2); 
        S(i,j,1,2)=-s1(1); 
        % 
        s2=[x(i+1,j+1)-x(i+1,j),y(i+1,j+1)-y(i+1,j)]; 
        S(i,j,2,1)=s2(2); 
        S(i,j,2,2)=-s2(1); 
        % 
        s3=[x(i,j+1)-x(i+1,j+1),y(i,j+1)-y(i+1,j+1)]; 
        S(i,j,3,1)=s3(2); 
        S(i,j,3,2)=-s3(1); 
        % 
        s4=[x(i,j)-x(i,j+1),y(i,j)-y(i,j+1)]; 
        S(i,j,4,1)=s4(2); 
        S(i,j,4,2)=-s4(1);     
    end % of j loop 
end % of i loop 
end % fsidevectors 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [U,hL,hR,vxL,vxR]=finitialU(xcen,ycen) 
% This sets up a dam break problem in the x direction (uniform in y). 
% Verification: Verified.  
% Inserts initial U values in computational cells and inserts ghost cells 
% all around with zero entries so U is 3x(NI2)x(NJ+2) so computational 
% indices go from i=j=1 to i=NI+1, j=NJ+1. 
[m,n]=size(xcen); % xcen has ghost cells around it 
NI=m-2; % number of computational cells in x direction 
NJ=n-2; % number of computational cells in y direction 
% 
d=3;  % number of PDEs 
U=zeros(d,NI+2,NJ+2); % create correct sized array with ghost cells 
xmidpoint=(xcen(2,2)+xcen(NI+1,2))/2; % middle x value of computational domain 
g=9.81; 
% 
% 1D Dam break.  Left and right states are either side  
% of the centre of computational domain.  Values are constant in y. 
% Left states 
hL=10.0; % left water depth [m] 
% disp('for verification') 
% hL=2/g; 
vL=[0,0]; % left water velocity vector [m/s] 
vxL=vL(1); 
phiL=hL*g; % left geopotential 
% Right states 
hR=5.0; % right water depth [m] 
% disp('for verification') 
% hR=1/g; 
vR=[0,0]; % right water velocity vector [m/s] 
vxR=vR(1); 
phiR=hR*g; % right geopotential 
% 
% Fill in U values for each computational cell 
for i=2:NI+1 
    for j=2:NJ+1 
        if (xcen(i,j)<=xmidpoint) 
           U(1,i,j)=phiL; 
           U(2,i,j)=phiL*vL(1); 
           U(3,i,j)=phiL*vL(2); 
        else 
           U(1,i,j)=phiR; 
           U(2,i,j)=phiR*vR(1); 
           U(3,i,j)=phiR*vR(2); 
        end % of if 
    end % of j loop 
end % of i loop 
% 
end % finitialU 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function U=fbcs(U) 
% Implements boundary conditions using ghost cells index by i=1, i=NI+2 
% j=1, j=NJ+2. 
[d,m,n]=size(U); 
NI=m-2; 
NJ=n-2; 
% 
% Zero gradient everywhere; 
U(:,1,:)=U(:,2,:); % left 
U(:,NI+2,:)=U(:,NI+1,:); % right 
U(:,:,1)=U(:,:,2); % bottom 
U(:,:,NJ+2)=U(:,:,NJ+1); % top 
end % fbcs 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function Unext=fsolver(U,A,dt,S,Dinterfaces,Dcentres) 
% Implements the finite volume scheme for a time step 
[d,m,n]=size(U); % includes ghost cells 
NI=m-2; % number of computational cells in i direction 
NJ=n-2; % number of computational cells in j direction 
for i=2:NI+1; % loop over computational cells 
    for j=2:NJ+1; % loop over computational cells 
     % put the 4 side vectors for cell (i,j) into 1x2 arrays 
     SRight=[S(i,j,2,1),S(i,j,2,2)]; 
     SLeft=[S(i,j,4,1),S(i,j,4,2)];  
     SUp=[S(i,j,3,1),S(i,j,3,2)]; 
     SDown=[S(i,j,1,1),S(i,j,1,2)];  
     % Find the 4 interface fluxes for cell i,j. 
     %disp('finterfaceflux') 
     [IHRight,IHLeft,IHUp,IHDown]=finterfaceflux(U,i,j,Dinterfaces,Dcentres); 
     % 
     % Update U to Unext 
     for k=1:d 
     % compute total flux out of cell i for equation k 
       IHR=[IHRight(k,1),IHRight(k,2)]; % right flux vector 
       IHL=[IHLeft(k,1),IHLeft(k,2)];   % left flux vector 
       IHU=[IHUp(k,1),IHUp(k,2)];       % upper flux vector 
       IHD=[IHDown(k,1),IHDown(k,2)];   % lower flux vector 
       %  
       totalfluxout=dot(IHR,SRight)+dot(IHL,SLeft)+dot(IHU,SUp)+dot(IHD,SDown); 
       %totalfluxout=dot(IHR,SRight)+dot(IHL,SLeft); %same results by just 
       % considering i direction fluxes which is correct for this problem 
       % Lax-Friedrichs scheme: 
       average=(U(k,i-1,j)+U(k,i+1,j)+U(k,i,j-1)+U(k,i,j+1))/4; 
       Unext(k,i,j)=average-(dt/A(i,j))*totalfluxout;  
     end % of k loop 
    end % of j loop 
 end % of i loop 
 end % fsolver 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [IHRight,IHLeft,IHUp,IHDown]=finterfaceflux(U,i,j,Dinterfaces,Dcentres) 
% Interface H values are based on linear extrapolation of neighbouring cell  
% centre data.  
% These may be found either by extrapolating U values or H components.  
% Since the PDEs are non-linear these approaches are not the same. 
% The latter appoach is analogous to the Lax-Friedrichs finite  
% difference method and is the better method as the former method produces 
% unwanted oscillations. 
% 
% Interface fluxes in each cell are denoted IHRight, IHLeft, IHUp, IHDown. 
[d,m,n]=size(U);   % d=number of rows of U (= number of PDEs = 3) 
NI=m-2;  % number of computational cells in i (x) direction 
NJ=n-2;  % number of computational cells in j (y) direction. 
IHRight=zeros(d,2); % array to store the right interface flux for a cell.  
                   % IHRight(,1) is the i component of H and 
                   % IHRight(,2) is the j component of H which is zero. 
IHLeft=zeros(d,2); 
IHUp=zeros(d,2); 
IHDown=zeros(d,2); 
% 
% First find H at cell centres. 
Hij=fgetHfromU(U,i,j); % Cell i,j 
Himinus1j=fgetHfromU(U,i-1,j); % Cell i-1,j 
Hiplus1j=fgetHfromU(U,i+1,j); % Cell i+1,j 
Hijminus1=fgetHfromU(U,i,j-1); % Cell i,j-1 
Hijplus1=fgetHfromU(U,i,j+1); % Cell i,j+1 
% 
% Side 1: Lower interface (i,j-1/2).  
GradH=(Hij-Hijminus1)/Dcentres(i,j,1); % gradients of H components 
IHDown=Hij-Dinterfaces(i,j,1)*GradH;  
% Side 2: Right interface (i+1/2,j).  
GradH=(Hiplus1j-Hij)/Dcentres(i,j,2); % gradients of H components 
IHRight=Hij+Dinterfaces(i,j,2)*GradH;  
% 
% Side 3: Upper interface (i,j+1/2).  
GradH=(Hijplus1-Hij)/Dcentres(i,j,3);  
IHUp=Hij+Dinterfaces(i,j,3)*GradH; 
% Side 4: Left interface (i-1/2,j).  
GradH=(Hij-Himinus1j)/Dcentres(i,j,4);  
IHLeft=Hij-Dinterfaces(i,j,4)*GradH; 
% 
end % finterfaceflux 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function H=fgetHfromU(U,i,j) 
% U = [   phi 
%       phi*vx 
%       phi*vy ] 
% 
% H=<     phi*vx    ,  phi*vy             
%    phi*vx^2+phi^2/2,  phi*vx*vy                    
%         phi*vx*vy ,  phi*vy^2+phi^2/2 > 
% 
[d,m,n]=size(U); 
H=zeros(d,2); % array holding flux density vector components in cell i,j 
% 
% The following could easily be made more efficient but we write out in  
% full for clarity. 
phi=U(1,i,j); 
vx=U(2,i,j)/phi; 
vy=U(3,i,j)/phi; 
% 
H(1,1)=phi*vx; 
H(1,2)=phi*vy; 
% 
H(2,1)=phi*vx*vx+phi*phi/2; 
H(2,2)=phi*vx*vy; 
% 
H(3,1)=phi*vx*vy; 
H(3,2)=phi*vy*vy+phi*phi/2; 
end % fgetHfromU 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function [Dinterfaces,Dcentres]=fcelldistances(x,y,xcen,ycen)  
% Dinterfaces holds the distances from centre of cell i,j to its 4 interfaces 
% measured along the lines joining the centres of neighbouring cells. 
% Note that if a lines are vertical then a special treatment is needed as 
% gradients are infinite (this will be the case with a Cartesian mesh). 
% Verification:  NOT VERIFIED 
% Dcentres holds the distsnces from the centre of cell i, j to the centres 
% of its 4 neighbouring cells. 
% 
[M,N]=size(xcen); % includes ghost cells 
NI=M-2; % number of computaational cells in i direction 
NJ=N-2; % number of computaational cells in j direction 
Dinterfaces=zeros(size(xcen)); % includes ghost cells 
Dcentres=zeros(size(xcen)); % includes ghost cells 
% 
tol=0.0001; % tolerance for special treatment (vertical lines) 
%  
for i=2:NI+1 % loop over computational cells 
    for j=2:NJ+1 % loop over computational cells 
        % Side 1: 
        a=xcen(i,j-1); % x coord of centre of cell (i,j-1) 
        b=ycen(i,j-1); % y coord of centre of cell (i,j-1) 
        c=xcen(i,j); % x coord of centre of cell (i,j) 
        d=ycen(i,j); % y coord of centre of cell (i,j) 
        e=x(i+1,j); % x coord of lower right corner of cell (i,j) 
        f=y(i+1,j); % y coord of lower right corner of cell (i,j) 
        g=x(i,j); % x coord of lower left corner of cell (i,j) 
        h=y(i,j); % y coord of lower left corner of cell (i,j) 
        % 
        % check lines are not vertical 
        if (abs(c-a) >= tol && abs(g-e) >= tol) % neither line is vertical 
           g1=(d-b)/(c-a); % gradient of line joining cell centres 
           g2=(h-f)/(g-e); % gradient of intersected cell side line 
           xi=(f-b+a*g1-e*g2)/(g1-g2); yi=b+g1*(xi-a); % intersection point     
        elseif (abs(c-a) < tol) % line joining cell centres is vertical 
          g2=(h-f)/(g-e); % gradient of intersected cell side line 
          xi=a; yi=f+g2*(xi-e); % intersection point     
        else % line along cell side is vertical 
          g1=(d-b)/(c-a); % gradient of line joining cell centres 
          xi=e; yi=b+g1*(xi-a); % intersection point               
          end % of if   
        Dinterfaces(i,j,1)=sqrt((xi-a)^2+(yi-b)^2); 
        %  
        % vector between centres of cells (i,j) and (i,j-1). 
        vec=[xcen(i,j)-xcen(i,j-1),ycen(i,j)-ycen(i,j-1)];                                               
        Dcentres(i,j,1)=norm(vec); % distance 
        % 
        % Side 2: 
        a=xcen(i,j); % x coord of centre of cell (i,j) 
        b=ycen(i,j); % y coord of centre of cell (i,j) 
        c=xcen(i+1,j); % x coord of centre of cell (i+1,j) 
        d=ycen(i+1,j); % y coord of centre of cell (i+1,j) 
        e=x(i+1,j); % x coord of lower right corner of cell (i,j) 
        f=y(i+1,j); % y coord of lower right corner of cell (i,j) 
        g=x(i+1,j+1); % x coord of upper right corner of cell (i,j) 
        h=y(i+1,j+1); % y coord of upper right corner of cell (i,j) 
        % 
        % check lines are not vertical 
        if (abs(c-a) >= tol && abs(g-e) >= tol) % neither line is vertical 
           g1=(d-b)/(c-a); % gradient of line joining cell centres 
           g2=(h-f)/(g-e); % gradient of intersected cell side line 
           xi=(f-b+a*g1-e*g2)/(g1-g2); yi=b+g1*(xi-a); % intersection point     
        elseif (abs(c-a) < tol) % line joining cell centres is vertical 
          g2=(h-f)/(g-e); % gradient of intersected cell side line 
          xi=a; yi=f+g2*(xi-e); % intersection point     
        else % line along cell side is vertical 
          g1=(d-b)/(c-a); % gradient of line joining cell centres 
          xi=e; yi=b+g1*(xi-a); % intersection point               
        end % of if   
        Dinterfaces(i,j,2)=sqrt((xi-a)^2+(yi-b)^2); 
        % 
        % vector between centres of cells (i,j) and (i+1,j). 
        vec=[xcen(i+1,j)-xcen(i,j),ycen(i+1,j)-ycen(i,j)];                                                  
        Dcentres(i,j,2)=norm(vec); % distance 
        % 
        % Side 3: 
        a=xcen(i,j); % x coord of centre of cell (i,j) 
        b=ycen(i,j); % y coord of centre of cell (i,j) 
        c=xcen(i,j+1); % x coord of centre of cell (i,j+1) 
        d=ycen(i,j+1); % y coord of centre of cell (i,j+1) 
        e=x(i+1,j+1); % x coord of upper right corner of cell (i,j) 
        f=y(i+1,j+1); % y coord of upper right corner of cell (i,j) 
        g=x(i,j+1); % x coord of upper left corner of cell (i,j) 
        h=y(i,j+1); % y coord of upper left corner of cell (i,j) 
        % 
        % check lines are not vertical 
        if (abs(c-a) >= tol && abs(g-e) >= tol) % neither line is vertical 
           g1=(d-b)/(c-a); % gradient of line joining cell centres 
           g2=(h-f)/(g-e); % gradient of intersected cell side line 
           xi=(f-b+a*g1-e*g2)/(g1-g2); yi=b+g1*(xi-a); % intersection point     
        elseif (abs(c-a) < tol) % line joining cell centres is vertical 
          g2=(h-f)/(g-e); % gradient of intersected cell side line 
          xi=a; yi=f+g2*(xi-e); % intersection point     
        else % line along cell side is vertical 
          g1=(d-b)/(c-a); % gradient of line joining cell centres 
          xi=e; yi=b+g1*(xi-a); % intersection point               
        end % of if   
        Dinterfaces(i,j,3)=sqrt((xi-a)^2+(yi-b)^2); 
        % 
        % vector between centres of cells (i,j) and (i,j+1). 
        vec=[xcen(i,j+1)-xcen(i,j),ycen(i,j+1)-ycen(i,j)];                                                  
        Dcentres(i,j,3)=norm(vec); % distance 
        % 
        % Side 4: 
        a=xcen(i-1,j); % x coord of centre of cell (i-1,j) 
        b=ycen(i-1,j); % y coord of centre of cell (i-1,j) 
        c=xcen(i,j); % x coord of centre of cell (i,j) 
        d=ycen(i,j); % y coord of centre of cell (i,j) 
        e=x(i,j); % x coord of lower left corner of cell (i,j) 
        f=y(i,j); % y coord of lower left corner of cell (i,j) 
        g=x(i,j+1); % x coord of upper left corner of cell (i,j) 
        h=y(i,j+1); % y coord of upper left corner of cell (i,j) 
        % 
        % check lines are not vertical 
        if (abs(c-a) >= tol && abs(g-e) >= tol) % neither line is vertical 
           g1=(d-b)/(c-a); % gradient of line joining cell centres 
           g2=(h-f)/(g-e); % gradient of intersected cell side line 
           xi=(f-b+a*g1-e*g2)/(g1-g2); yi=b+g1*(xi-a); % intersection point     
        elseif (abs(c-a) < tol) % line joining cell centres is vertical 
          g2=(h-f)/(g-e); % gradient of intersected cell side line 
          xi=a; yi=f+g2*(xi-e); % intersection point     
        else % line along cell side is vertical 
          g1=(d-b)/(c-a); % gradient of line joining cell centres 
          xi=e; yi=b+g1*(xi-a); % intersection point               
        end % of if   
        Dinterfaces(i,j,4)=sqrt((xi-a)^2+(yi-b)^2); 
        % 
        % vector between centres of cells (i-1,j) and (i,j). 
        vec=[xcen(i,j)-xcen(i-1,j),ycen(i,j)-ycen(i-1,j)];   
        Dcentres(i,j,4)=norm(vec); % distance 
    end % j loop 
end % i loop 
end % fcelldistances 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=fdisplay(in) 
 % displays matrices 'correctly' - WORKS 
 % Using index i for x which are columns and j for y which are rows in  
 % reverse order so this routine displays them in this fashion 
 [NX,NY]=size(in); 
 temp=transpose(in); 
 out=temp(NY:-1:1,:); 
end % fdisplay 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function dt=fcalcdt(A,S,U) 
% Verification:  needs checking, looks wrong 
% Finds allowable time step dt using heuristic formula. 
% 0<F<1 is a tunable safety factor. 
[m,n]=size(A);  % A and S include ghost cells all around. 
NI=m-2; % number of computational cells in i direction 
NJ=n-2; % number of computational cells in j direction 
dt=9999999;  % large value of dt to start 
F=0.5; 
 for i=2:NI+1 
    for j=2:NJ+1 
        % put side vectors into 1x2 arrays   
        SRight=[S(i,j,2,1),S(i,j,2,2)]; 
        SUp=[S(i,j,3,1),S(i,j,3,2)]; 
        % 
        phi=U(1,i,j); 
        cel=sqrt(phi); 
        vx=U(2,i,j)/phi; 
        vy=U(3,i,j)/phi; 
        v=[vx,vy]; % velocity vector 
        % 
        vR=dot(v,SRight); % velocity component through right cell side 
        vU=dot(v,SUp); % velocity component through upper cell side 
        b=cel*norm(SRight); 
        c=cel*norm(SUp); 
        speedR=max(abs(vR+b),abs(vR-b));   % CHECK 
        speedU=max(abs(vU+c),abs(vU-c));    
        dum=A(i,j)/max(speedR,speedU); 
        % take smallest value of dt 
        if (dum<dt) 
           dt=dum; 
        end 
    end % j loop 
 end % i loop 
dt=F*dt; % multiply dt by safety factor F 
end  % fcalcdt 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function fdrawmesh(x,y) 
% Description:  Draws a structured 2D finite volume mesh and fills in any 
% solid cells (not coded here). 
% Date structures: 
% The computational mesh has NI cells in the i direction and  
% NJ cells in the j direction.  x and y are surrounded by ghost cells so 
% that x and y are both (NI+3)x(NJ+3). 
% The x and y coordinates of the lower left hand corner of cell (i,j) are 
% held in arrays x and y respectively.  In 
% this way the 4 vertices of cell (i,j) are (x(i),y(j)), (x(i+1),y(j)), 
% (x(i+1),y(j+1)) and (x(i),y(j+1)).  solid is an NI+2 by NJ+2 array which  
% flags solid cells.  If cell (i,j) is solid then solid(i,j)=1 otherwise 
% solid(i,j)=0. 
%  
[m,n]=size(x); 
NI=m-3; % number of computational cell in i direction 
NJ=n-3; % number of computational cell in j direction 
% 
% Plot the mesh 
hold on   % keeps all plots on the same axes 
% draw lines in the i direction 
for i=2:m-1 
    plot(x(i,:),y(i,:)) 
end 
% draw lines in the j direction 
for j=2:n-1 
    plot(x(:,j),y(:,j))    
end 
title('computational mesh') 
xlabel('x') 
ylabel('y') 
% % Fill in solid cells 
% for i=1:NI+1 
%     for j=2:NJ+1 
%         if (solid(i,j)==1) 
%             solidx=[x(i,j),x(i+1,j),x(i+1,j+1),x(i,j+1),x(i,j)]; 
%             solidy=[y(i,j),y(i+1,j),y(i+1,j+1),y(i,j+1),y(i,j)]; 
%             fill(solidx,solidy,'k')       
%         end 
%     end % of j loop 
% end % of i loop 
end % fdrawmesh 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function fplotresults(xcen,ycen,U) 
% Ghost cells have been removed. 
% Plots results on a structured mesh. 
% 
[NI,NJ]=size(xcen); 
midj=floor((NJ+1)/2);  % middle j index 
g=9.81;  % acceleration due to gravity 
% plot in x direction at middle j index 
  subplot(3,1,1), plot(xcen(:,midj),U(1,:,midj)/g) 
  title('numerical results: plot of water depth at middle j index') 
  xlabel('x [m]') 
  ylabel('h [m]') 
  % 
  subplot(3,1,2), plot(xcen(:,midj),U(2,:,midj)./U(1,:,midj)) 
  title('numerical results: plot of water speed in x at middle j index') 
  xlabel('x [m]') 
  ylabel('vx [m/s]') 
  % 
  subplot(3,1,3), plot(xcen(:,midj),U(3,:,midj)./U(1,:,midj)) 
  title('numerical results: plot of water speed in y at middle j index') 
  xlabel('x [m]') 
  ylabel('vy [m/s]') 
  
% more general routines: 
%   subplot(3,1,1), contour(xcen,ycen,U(1,:,:)/g) 
%   title('numerical results: contour plot of water depth') 
%   xlabel('x [m]') 
%   ylabel('y [m]') 
%   % 
%   subplot(3,1,2), contour(xcen,ycen,U(2,:,:)/U(1,:,:)) 
%   title('numerical results: contour plot of water speed in x') 
%   xlabel('x [m]') 
%   ylabel('y [m]') 
%   % 
%   subplot(3,1,3), contour(xcen,ycen,U(3,:,:)/U(1,:,:)) 
%   title('numerical results: contour plot of water speed in y') 
%   xlabel('x [m]') 
%   ylabel('y [m]') 
% 
end % fplotresults 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% The subfunctions below are for finding the exact solution to the 1D dam 
% break problem and may be commented out if not required 
function [eh,ev]=fexactSWE1D(runtime,xcen,hL,hR,uL,uR) 
% Author: Clive Mingham 
% Date: Nov 24 2011 
% Description: Gives exact solution for water depth and water speed to the  
% dam break (Riemann) problem for the 1D shallow water equations. 
%  
% Left and right states defined are either side of x = 0.0  
% 
% Verification:  Looks right for type 2. 
% 
g=9.81; 
time=runtime; 
% 
[NI,NJ]=size(xcen); 
M=floor((NJ+1)/2); % middle j value 
% xRiemann is true x coordinate of the initial Riemann problem 
xRiemann=(xcen(1,M)+xcen(NI,M))/2; % middle x coordinate 
% 
% Initial states either side of x = 0 so results must be translated to  
% fit the problem domain. 
% hL is left height [m] 
% hR is right height [m] 
% uL is left speed [m/s] 
% uR is right speed [m/s] 
% 
phiL=g*hL; % left geopotential 
phiR=g*hR; % right geopotential 
aL=sqrt(phiL); % left celerity 
aR=sqrt(phiR); % right celerity 
% 
% Initial large values of wave fronts 
x1=9999; 
x2=9999; 
x3=9999; 
x4=9999; 
% The non-linear equation for phis (star region) is solved using the   
% Newton-Raphson iteration method, an initial guess is: 
phiold=(0.5*(sqrt(phiL)+sqrt(phiR))+0.25*(uL-uR))^2; 
phis=fNR(phiold,phiL,phiR,uL,uR); % finds phi in the star region 
as=sqrt(phis); 
% us is the particle speed in the star region (constant) 
us=0.5*(uL+uR)+0.5*(ffr(phis,phiR)-ffl(phis,phiL)); 
% Having found phis we know which way shocks and expansion waves are 
% moving.  Then their speed can be found.  sL, sR are the shock speeds. 
% 
if (phis<=phiL && phis<=phiR) 
   %disp('type1 left and right expansion waves') 
   type=1; 
else    
  if (phis<=phiL && phis>phiR) 
     %disp(' type2 left expansion wave and right shock wave') 
     mr=sqrt(0.5*phiR*phis*(phiR+phis));  
     sR=uR+mr/phiR;            
     type=2;      
  else  
    if (phis>phiL && phis<=phiR) 
       %disp(' type3 left shock wave and right expansion wave') 
       ml=sqrt(0.5*phiL*phis*(phiL+phis)); 
       sL=uL-ml/phiL; 
       type=3; 
    else 
      %disp('type4 left shock wave and right shock wave') 
      ml=sqrt(0.5*phiL*phis*(phiL+phis)); 
      sL=uL-ml/phiL; 
      mr=sqrt(0.5*phiR*phis*(phiR+phis));  
      sR=uR+mr/phiR; 
      type=4; 
    end 
  end 
end 
% Find positions x1, x2 .. of the wave fronts for the different cases 
if (type==1) 
   x1=(uL-aL)*time; 
   x2=(us-as)*time; 
   x3=(us+as)*time; 
   x4=(uR+aR)*time; 
   right=x4; 
else     
  if (type==2) 
     x1=(uL-aL)*time; 
     x2=(us-as)*time; 
     x3=sR*time; 
     right=x3;      
  else 
    if (type==3) 
       x1=sL*time;  
       x2=(us+as)*time; 
       x3=(uR+aR)*time; 
       right=x3; 
    else    
      x1=sL*time; 
      x2=sR*time; 
      right=x2; 
    end 
  end 
end 
% Generate exact solutions 
xx=xcen(:,M); % x coords of cell centres on middle row 
ev=zeros(NI,1); % exact speed  
eh=zeros(NI,1); % exact depth 
% 
for i=1:NI 
    x=xx(i)-xRiemann; % x coord using xRiemann as origin 
    if (type==1) 
       if (x<x1)       
          ev(i)=uL; 
          ephi(i)=phiL; 
       else    
         if (x<x2) 
            ev(i)=fleu(x,time,aL,uL);  
            ephi(i)=flef(x,time,aL,uL); 
         else 
           if (x<x3) 
              ev(i)=us;  
              ephi(i)=phis; 
           else 
             if (x<x4) 
                ev(i)=freu(x,time,aR,uR);  
                ephi(i)=fref(x,time,aR,uR); 
             else 
               ev(i)=uR;  
               ephi(i)=phiR; 
             end      
           end 
         end 
       end 
     else 
       if (type==2) 
          %disp('calculating type 2 results') 
          if (x<x1) 
             ev(i)=uL;  
             ephi(i)=phiL; 
          else       
            if (x<x2) 
               ev(i)=fleu(x,time,aL,uL);  
               ephi(i)=flef(x,time,aL,uL); 
            else 
              if (x<x3) 
                 ev(i)=us; 
                 ephi(i)=phis; 
              else 
                ev(i)=uR; 
                ephi(i)=phiR; 
              end 
            end       
          end    
        else 
          if (type==3) 
             if (x<x1) 
                ev(i)=uL;  
                ephi(i)=phiL; 
             else       
               if (x<x2) 
                  ev(i)=us;  
                  ephi(i)=phis; 
               else 
                 if (x<x3) 
                    ev(i)=freu(x,time,aR,uR);  
                    ephi(i)=fref(x,time,aR,uR); 
                 else 
                   ev(i)=uR;  
                   ephi(i)=phiR; 
                 end 
               end 
             end    
          else     
            if (x<x1) 
               ev(i)=uL;  
               ephi(i)=phiL; 
            else       
              if (x<x2) 
                 ev(i)=us; 
                 ephi(i)=phis; 
              else 
                ev(i)=uR;  
                ephi(i)=phiR; 
              end 
            end 
          end    
       end 
    end 
end % of i loop 
eh=ephi/g; 
% 
% disp('see plots of exact results') 
% subplot(2,1,1), plot(xx,eh) 
% title('Exact solution to Riemann problem for 1D SWE') 
% xlabel('x [m]') 
% ylabel('water depth [m]') 
% subplot(2,1,2), plot(xx,ev) 
% title('Exact solution to Riemann problem for 1D SWE') 
% xlabel('x [m]') 
% ylabel('water speed [m/s]') 
end % of function fexactSWE1D 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=fNR(x,phiL,phiR,uL,uR) 
% Newton-Raphson method for solving a non-linear equation f(x)=0 
% initially x is the initial guess. 
tol=0.00001;  
error=0.0001; 
itermax=100; 
itertot=0; 
% Start iteration  
for j=1:itermax 
         itertot=itertot+1; 
         if (j==itermax) 
            disp('WARNING: max iterations achieved without convergence') 
         end 
         if (abs(fdfunc(x,phiL,phiR))<tol) 
             disp('DERIVATIVE TOO CLOSE TO ZERO, NR METHOD FAILS') 
             break 
         end 
         xnew=x-ffunc(x,phiL,phiR,uL,uR)/fdfunc(x,phiL,phiR); 
         reldif=abs((xnew-x)/x); 
         if (reldif<error) 
            disp('iteration finished') 
            out=xnew; 
            break 
         else 
            x=xnew; 
         end                
end % j loop 
end % function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=ffunc(x,phiL,phiR,uL,uR) 
% function is f(x) = 0 
  out=ffl(x,phiL)+ffr(x,phiR)+uR-uL; 
end % function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=fdfunc(x,phiL,phiR) 
% f'(x) in f(x) = 0 
  out=fdfl(x,phiL)+fdfr(x,phiR); 
end % function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=ffl(x,phiL) 
  if (x<=phiL) 
     out=2*(sqrt(x)-sqrt(phiL)); 
  else 
     out=(x-phiL)*sqrt(0.5*(x+phiL)/(x*phiL)); 
  end 
end % function 
%%%%%%%%%%%%%%%%%%%%%**%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=ffr(x,phiR)  
  if (x<=phiR) 
     out=2*(sqrt(x)-sqrt(phiR)); 
  else 
    out=(x-phiR)*sqrt(0.5*(x+phiR)/(x*phiR)); 
  end 
end % function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=fdfl(x,phiL) 
  dL=sqrt(0.5*(x+phiL)/(x*phiL)); 
  if (x<=phiL) 
     out=x^(-0.5); 
  else 
    out=dL-0.25*(x-phiL)/(dL*x^2); 
  end 
end % function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=fdfr(x,phiR) 
  dr=sqrt(0.5*(x+phiR)/(x*phiR)); 
  if (x<=phiR) 
     out=x^(-0.5); 
  else 
    out=dr-0.25*(x-phiR)/(dr*x^2); 
  end 
end % function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=fleu(x,t,aL,uL) 
% u value in left expansion wave      
  out=(uL+2*aL+2*x/t)/3; 
end % function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=freu(x,t,aR,uR) 
% u value in right expansion wave 
  out=(uR-2*aR+2*x/t)/3; 
end % function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=flef(x,t,aL,uL) 
% phi in left expansion wave 
  out=((uL+2*aL-x/t)/3)^2; 
end % function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function out=fref(x,t,aR,uR) 
% phi in right expansion wave 
  out=((-uR+2*aR+x/t)/3)^2; 
end % function 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
 

